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ABSTRACT 

It has been suggested that the transition between magnetorotationally active and dead zones 
in protoplanetary disks should be prone to the excitation of vortices via Rossby wave instability 
(RWI). However, the only numerical evidence for this has come from alpha disk models, where 
the magnetic field evolution is not followed, and the effect of turbulence is parametrized by Lapla- 
cian viscosity. We aim to establish the phenomenology of the flow in the transition in 3D resistive- 
magnetohydrodynamical models. We model the transition by a sharp jump in resistivity, as expected 
in the inner dead zone boundary, using the PENCIL CODE to simulate the flow. We find that vortices 
are readily excited in the dead side of the transition. We measure the mass accretion rate finding 
similar levels of Reynolds stress at the dead and active zones, at the a. « 10~ 2 level. The vortex sits 
in a pressure maximum and does not migrate, surviving until the end of the simulation. A pres- 
sure maximum in the active zone also triggers the RWI. The magnetized vortex that results should 
be disrupted by parasitical magneto-elliptic instabilities, yet it subsists in high resolution. This sug- 
gests that either the parasitic modes are still numerically damped, or that the RWI supplies vorticity 
faster than they can destroy it. We conclude that the resistive transition between the active and dead 
zones in the inner regions of protoplanetary disks, if sharp enough, can indeed excite vortices via 
RWI. Our results lend credence to previous works that relied on the alpha-disk approximation, and 
caution against the use of overly reduced azimuthal coverage on modeling this transition. 



1. INTRODUCTION 

The formation of planets remains one of the most chal- 
lenging problems of contemporary astrophysics. The 
current paradigm in planet formation theory describes 
a hierarchical growth of solid bodies, from interstel- 
lar dust grains to rocky planetary cores (Safronov 1969; 
Lyttleton 1972; Goldreich & Ward 1973; Youdin & Shu 
2002). A particularly difficult phase in the process is the 
growth from centimeter sized pebbles and meter-sized 
boulders to planetary embryos the size of our Moon or 
Mars. Objects in the pebble to boulder range are ex- 
pected to drift inward extremely rapidly in a protoplan- 
etary disk, so that they would generally fall into the cen- 
tral star well before larger bodies can form by simple 
accumulation (Weidenschilling 1977; Brauer et al. 2008). 

Ways to bypass this problem have focused on inho- 
mogeneities in the flow, in order to trap particles in 
their migrating path. Cuzzi et al. (2008) proposed a 
model in which mm-sized particles are trapped in the 
smallest eddies in the flow, forming "sandpile" plan- 
etesimals in the 10-100 km range (though that model 
has been criticized by Chang & Oishi 2010 and Pan et al. 
2011). Particles may also be trapped in mesoscale "zonal 
flows" (Lyra et al. 2008a; Johansen et al. 2009; Simon et 
al. 2012) that are local inversions in the angular velocity 
profile, brought about by spatial variations in magnetic 
pressure. The particles themselves can give rise to the 
necessary inhomogeneities, as their migrating stream- 
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ing flow develops into a traffic-jam instability (Youdin & 
Shu 2002; Youdin & Goodman 2005; Youdin & Johansen 
2007), leading to intense particle clumping (Johansen & 
Youdin 2007) and subsequent planetesimal formation at 
the Ceres-mass range (Johansen et al. 2007). It has also 
been recently proposed that icy planetesimals may form 
from direct coagulation (Okuzumi et al. 2012) due to the 
enhanced sticking properties of ices. 

Another process has been suggested, that combines 
several of the advantages (as well as many of the prob- 
lems) of the scenarios described above, and is the sub- 
ject of this work. Turbulence in the largest scales of the 
flow, in the form of large scale vortices, has been inde- 
pendently proposed by Barge & Sommeria (1995) and 
Tanga et al. (1996) as fast routes for planet formation, 
for two main reasons. First, vortices are equilibrium 
solutions of the Navier-Stokes equations, and thus are 
persistent structures in hydrodynamic flows, as seen in 
the Great Red Spot of Jupiter, a remarkable high pres- 
sure vortex stable since first spotted, over three hun- 
dred years ago (Hooke 1665; Cassini 1666) 4 . The second 
is that the equilibrium is geostrophic, i.e, between the 
Coriolis force and the pressure gradient force. As solids 
do not feel the pressure force, the Coriolis force will lead 
them out of the vortex if cyclonic and into the eye if an- 
ticyclonic. As the shear enforces that only anti cyclonic 
vortices persist (Marcus 1993; Adams & Watkins 1995; 
Bracco et al. 1999; Godon & Livio 1999), this becomes a 
very effective mechanism to concentrate solid particles 
(Klahr 2006), as also observed in numerical simulations 
(Godon & Livio 2000; Johansen et al. 2004; Fromang & 

4 The reader is refereed to the fascinating account of Falorni (1987) 
on the history of the discovery of the Spot, where the author debates 
the claims of primacy to Hooke or Cassini. 
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Nelson 2005; Inaba & Barge 2006). It was further shown 
by Lyra et al. (2008b, 2009) that in the limit where turbu- 
lence in the vortex core is absent and the fragmentation 
of particles is ignored, the concentration of solids easily 
reaches the conditions necessary to gravitationally col- 
lapse them into planets. 

Nevertheless, as exciting as the vortex hypothesis 
may be, for a long time no plausible mechanism for 
the formation and sustenance of such storm systems 
in disks could be found. It is well known that three- 
dimensional vortices fall prey to elliptical instabilities 
(Bayly 1986; Pierrehumbert 1986; Kerswell 2002; Lesur 
& Papaloizou 2009), a general name given to parasitic 
instabilities of closed elliptical streamlines (see e.g. the 
review of Kerswell 2002). Fluids in rigid rotation sup- 
port a spectrum of stable inertial waves, the simplest 
ones being circularly polarized transverse plane waves 
oscillating at twice the base frequency (see e.g. Chan- 
drasekhar 1961). Strain is introduced when the mo- 
tion passes from circular to elliptical, and some three- 
dimensional modes find resonance with the underlying 
strain field. Because of its intrinsic three-dimensional 
character, the absence of elliptic instability modifies the 
turbulent energy cascade from direct to inverse in 2D 
(Batchelor 1969). The result is that the series of vis- 
cous mergings that turn small and mesoscale eddies 
into large vortices in the integral scale is not expected 
to occur in systems that significantly depart from two- 
dimensionality. In other words, vortices do not sponta- 
neously develop in three dimensions from initial small 
scale noise as they do in two dimensions, and must 
therefore be the result of some instability. Because the 
Keplerian shear provides stabilization of linear axisym- 
metric disturbances at all Reynolds numbers (Rayleigh 
criterion; Strutt 1880, 1916), in non-magnetized disks 
any such instability must be either non-axisymmetric, 
non-linear, or rely on a modification of the angular ve- 
locity profile to circumvent the stabilizing effect of the 
shear. 

Based on these ideas, two main processes have been 
proposed for the formation of vortices in disks. One is 
the non-linear radial convection process that has been 
referred to as global baroclinic instability (Klahr & Bo- 
denheimer 2003; Klahr 2004), subcritical baroclinic in- 
stability (Lesur & Papaloizou 2009; Paardekooper et al. 

2010) , and simply baroclinic instability (Lyra & Klahr 

2011) . This shall not be dealt with in this barotropic pa- 
per. We will focus here on what has become known as 
Rossby wave instability (RWI; Lovelace et al. 1999; Li et 
al. 2000, 2001; Umurhan 2010). 

The RWI relies on a modification of the angular ve- 
locity profile, brought about by a local extremum of an 
entropy-modified potential vorticity quantity (or simply 
the potential vorticity in the case of a barotropic flow). 
The extremum in potential vorticity launches inertial- 
acoustic waves, akin to Rossby waves in planetary at- 
mospheres, and traps modes in its co-rotational singu- 
larity (see fig 1 of Meheut et al. 2010). The mechanism 
was first discussed by Lovelace & Hohlfeld (1978) in 
the context of self-gravitating galactic disks, and men- 
tioned en passant by Toomre (1981), who called them 
"edge modes". The instability was discussed again 
by Papaloizou & Pringle (1984, 1985) in the context 
or so-called "slender torus" models, hot disks around 



quasars modeled locally in radius. The first non-linear 
numerical simulation of the instability, done by Haw- 
ley (1987), found that the trapped non-axisymmetric 
modes evolved into anticyclonic vortices, but could not 
follow the calculation until full saturation. The inter- 
est on the instability waned after the re-discovery of 
the magneto-rotational instability (MRI; Velikhov 1959; 
Chandrasekhar 1960) by Balbus & Hawley (1991) and 
its rise to paradigmatic status as the source of angu- 
lar momentum transport and turbulence in disks. The 
RWI was re-analyzed by Lovelace et al. (1999), who 
expanded the linear analysis of Papaloizou & Pringle 
(1984) to non-barotropic disks, and gave the mechanism 
its modern name. Subsequent work showed that a local- 
ized bump in surface density or pressure would cause 
the instability (Li et al. 2000), and numerical simula- 
tions (Li et al. 2001) confirmed its saturation into vor- 
tices, with most unstable azimuthal wavenumbers m=4- 
5. 

Varniere & Tagger (2006) raised the possibility that 
such surface density jumps would naturally occur at the 
boundaries between the magnetized and unmagnetized 
regions of accretion disks. In this boundary, because the 
unmagnetized region constitutes a "dead zone" to the 
MRI (Gammie 1996), there exists a transition in turbu- 
lent viscosity. The viscous torque at the transition has a 
component proportional to the negative of the viscosity 
gradient, so material is accelerated outward in the inner 
dead zone boundary (negative viscous gradient) and in- 
ward in the outer dead zone boundary (positive viscous 
gradient). This modifies the potential vorticity profile at 
these transitions, triggering the RWI. Because this com- 
ponent of the viscous torque is also proportional to the 
shear rate, this process occurs at the shear timescale, not 
on the much longer viscous timescale. 

This two-dimensional scenario, modeling the turbu- 
lent transition as jumps in alpha-viscosity (Shakura 
& Sunyaev 1973) was used by Inaba & Barge (2006) 
to argue for accumulation of particles in that region, 
and by Lyra et al. (2008b, 2009) to demonstrate their 
collapse into planets. The test of this scenario in 
three-dimensional, magneto-hydrodynamical calcula- 
tions has yet to be performed. That is the goal of this 
paper. 

In particular, we focus on three questions. First, 
whether the RWI will drive angular momentum trans- 
port. Previous works (Fleming & Stone 2003; Oishi & 
Mac Low 2009) show that although turbulence in the 
active upper layers drives waves that propagate into 
the dead zone, the amount of stress generated is too 
low to drive significant angular momentum transport, 
due to the low inertia of the upper active layers. In 
this work we investigate the effect of a radial, not ver- 
tical, resistive transition. As in the midplane the den- 
sities are much higher, perhaps the degree of angular 
momentum transport changes significantly. Moreover, 
strong anti-cyclonic vortices provide an exciting possi- 
bility. Vorticity generates spiral density waves (Heine- 
mann & Papaloizou 2009a,b, 2012) that would have the 
right velocity correlation to transport angular momen- 
tum outwards. It has been shown in local models of 
unmagnetized disks that by radiation of waves, vor- 
tices can maintain a relatively high accretion rate, at 
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the level of ol = 10~ 3 (Lesur & Papaloizou 2009; Lyra & 
Klahr 2011). Varniere & Tagger (2006) showed proof- 
of-concept alpha-disk models that the dead zone can be 
"revived" by vortices excited in the active/ dead bound- 
ary, but that still has to be confirmed by MHD simula- 
tions. 

The second question pertains to the width of the tran- 
sition. Lyra et al. (2009; see also Regaly et al. 2012) show 
that the maximum width where one expects a pressure 
bump to trigger the RWI is 2H, where H is the pressure 
scale height. The outer edge of the dead zone under- 
goes a gradual transition, where the density falls mono- 
tonically, until the point that it gets so thin that ioniz- 
ing agents (stellar X-rays or external cosmic rays) can 
penetrate all the way to the midplane. As this transi- 
tion should be very smooth, extending through tens of 
AU, we can exclude that as a possible location for RWI 
excitation. The inner edge, on the other hand, houses 
a more exciting possibility, since there the transition 
in ionization is sharp, set by the collisional ionization 
of potassium, that occurs at ps 900 K (e.g. Umebayashi 
1983; Turner & Drake 2009). 

The other question is whether the magnetic field 
would in any way inhibit the growth of the RWI. 
Dzyurkevich et al. (2010) find that a density bump de- 
velops at the inner edge of the dead zone, but in the 
timespan of their simulation, no appreciable growth of 
non-axisymmetric modes is observed. The same nega- 
tive result is found in the models of Kato et al. (2009, 
2010, 2012) who model inhomogeneous MRI, mimick- 
ing the effects of a dead zone. This is curious, since 
non-axisymmetric modes trapped in the pressure bump 
should develop counter-rotation, as the fluid element 
ahead of the bump is accelerated, and the fluid element 
behind it is decelerated (Hawley 1987). These nega- 
tive results could be due to the limited azimuthal ex- 
tent modeled by both works (tt/4 by Dzyurkevich et al. 
2010, local box in the works of Kato et al.), of stratifica- 
tion in the case of Dzyurkevich et al. (2010), or inhibition 
by the magnetic field. 

Although vortices have been reported in MRI-active 
disks (Fromang & Nelson 2005), the magneto-elliptic in- 
stability (Mizerski & Bajer 2009; Mizerski & Lyra 2012) 
was shown by Lyra & Klahr (2011) to be powerful 
enough to disrupt vortices otherwise stable in non- 
magnetized environments. In this work we present 
models where the dead zone is modeled by a resistive 
jump in a magnetized disk, and show that the RWI is 
excited at the dead side of the transition, leading to a gi- 
ant vortex that survives until the end of the simulation. 
The paper is organized as follows. In Sect 2 we present 
the model equations, and in Sect 3 the initial conditions. 
The results are shown in Sect 4, followed by a discussion 
and conclusions in Sect 5. 



2. THE MODEL 

2.1. Dynamical equations 

We perform three-dimensional MHD simulations in 
the cylindrical approximation, i.e., neglecting the disk 
vertical stratification and switching off gravity in that 
direction. The equations solved are 
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Fig. 1. — Snapshots of density (left), residual vorticity (middle) and 
Alfven speed (right) at the midplane in selected times at 3, 5, and 50 
reference orbits. 



dp 
di 
du 
1ft 
OA 



■-(u-V)p-pV u, (1) 
:-(w V)w- ^Vp- V<£ + (2) 

■uxB-n}ioJ (3) 

V=P c l- ( 4 ) 

where p is the density, u the velocity, A is the mag- 
netic potential, B = V x A is the magnetic field, / = 



4 



Lyra & Mac Low 



1.4 
1.2 

Q_ 

0.8 
0.6 

1.02 
1.01 

O 

"€ 1 00 

s 

0.99 

0.98 
2.0 



1.5 



3 



1.0 



0.5 

0.0 
0.12 

0.10 

0.08 

5 0.06 

0.04 

0.02 
0.00 



7 


Density 


\ \ >' / 


r/r (J =j 








trr =i _ 


h-^-Hf - ■ 


1 ■ ■ 1 ■ 1 < 1 : 

Az. Velocity \ 







/ T 

/ill / \ 


Vorticity 








Mag. Field 



0.5 



1.0 



1.5 



2.C 



r/t\, 



FIG. 2. — Vertical and azimuthal averages of density, azimuthal ve- 
locity, vorticity, and magnetic field strength at 1, 2, and 3 orbits, i.e., 
during linear growth of the MRI. The build-up of magnetic stress 
at the resistive transition at r, t = 1.2 gives rise to an inversion in 
azimuthal velocity, with super-Keplerian motion inwards and sub- 
Keplerian motion outwards, generating a dip in potential vorticity. 
The opposite occurs in the artificial peak of magnetic pressure at 
r = 0.7 coded in the initial condition, leading to a peak in potential 
vorticity. Both locations fulfill the conditions for triggering the RWI. 

jiQ- V x B is the current density, and p is the pressure. 
The equation of state is locally isothermal. The grav- 
itational potential <P = -GM^/r where G is the gravi- 
tational potential, M* is the stellar mass, and r is the 
cylindrical radius. The resistivity is a radial function of 
position. We use a smooth step function 



n{r) 



2 



1 + tanh f - 



Ar 



(5) 



in order to mimic the effect of a dead zone. The resistiv- 
ity passes from r\§ to zero over a width A r centered at 
an arbitrarily chosen distance . 

We solve the equations with the PENCIL CODE 5 
5 The code, including improvements done for the present work, is 




FIG. 3. — Time evolution of the (p-z averaged (r-dependent) den- 
sity (left panel), vertical vorticity (middle panel) and magnetic energy 
(right panel) for the fiducial model. A density enhancement is seen 
at the dead side of the resistive transition at r = 1.2, matching radial 
locations of lower vorticity. Curiously, other vortices are excited in the 
active zone. One exists at /- « 0.7 from t ?s 3Tg to f ~ 20Tb, and another 
next to the inner boundary from f fa 40To to f ~ 70Tb, after which it 
decays and restarts at t fs 85To. 



which integrates the evolution equations with sixth or- 
der spatial derivatives, and a third order Runge-Kutta 
time integrator. Sixth-order hyper-dissipation terms are 
added to Eq. (l)-Eq. (3), to provide extra dissipation 
near the grid scale, explained in Lyra et al. (2008a). They 
are needed because the high order scheme of the Pencil 
Code has little overall numerical dissipation (McNally 
et al. 2012). 

2.2. Initial Conditions 

We model a three-dimensional disk on a uniformly 
spaced mesh in cylindrical coordinates (r,cp,z), ranging 
over r=[0.4,2.0]ro and z=[-0.1,0.1]ro, where r$ is a refer- 
ence radius. We run models with azimuthal coverage 
L^=7i and 



=27T. The fiducial model with L^=7r has 
resolution [N r ,N^,N 2 ]=[192,384,64]. 

The density and sound speed are set as radial power- 
laws 



P = Po 



-<Jf> 



-sO 



-<? T 



(6) 



with q p = 1.5 and q T = 1.0. The initial angular velocity 
profile is corrected by the thermal pressure gradient 



rp or 



(7) 



where Q = J?o iff ?'o)~' ? with q = 1.5 is the Keplerian an- 
gular velocity. 

The magnetic field is set as a net vertical field, with 
four MPJ wavelengths resolved in the vertical range. 
The constraint Amri = 2nv A fl~ l = L 2 /4 translates into 
a radially varying field 

publicly available under a GNU open source license and can be down- 
loaded at http: / /www.nordita.org/ software/ pencil-code 
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i.e., a field falling as a 9/4 power-law. We use units such 
that 

GM, = r = Po = p = 1 (9) 

and omit these factors hereafter. Thus, v A0 « 8 x 10" 3 in 
code units. The reference sound speed is set at c s o=0.1. 
The dimensionless plasma beta parameter f> — 2c 2 /i> 2 
then ranges from 50 in the inner disk to 1250 in the outer. 

Noise is added to the velocities at the 10" 4 level. In this 
configuration, the MRI grows and saturates quickly in 3 
local orbits. We apply the noise and the magnetic field 
only in the radial range r =[0.6,1.8] to avoid growth of 
the instability near the boundaries. We use reflective 
boundaries, with a buffer zone of width 0.2 at each ra- 
dial border, that drives the quantities to the initial con- 
dition on a dynamical timescale. 

In the presence of resistivity, the excitation of the MRI 
is controlled by the Elsasser number 



V 



(10) 



which is the magnetic Reynolds number with velocity 
equal to the Alfven speed, and A the relevant magnetic 
length scale. We set the reference resistivity t]q so that 
the Elsasser number of the largest wavelength present 
in the box (i.e., A = L z ) is unity at tq, thus quenching the 
MRI outward of that radius. This constraint translates 
into j/o = L z v A0 = 1.6 x 10~ 3 . We place the resistivity 
jump at r v =1.2, so that at that location A = 0.75 < 1. 

3. RESULTS 

3.1. Fiducial model 

The fiducial model ranges n in azimuth, with a resolu- 
tion of [N^N^Nzh [192,384,64], and a sharp resistivity 

transition at r n — 1.2 of width Ar = 10~ 2 . The MRI starts 
from the inner disk and quickly saturates in 3 local or- 
bits. We show in Fig. 1 the state of the disk for density 
(left panels), residual vertical vorticity (vertical vortic- 
ity minus the vertical Keplerian vorticity, middle pan- 
els) and Alfven speed (right panels) in selected snap- 
shots. The dashed line marks the boundary between 
MRI-active and dead zones. We see that the magnetic 
energy is well confined to the active zone, with only 
some very minor diffusion into the dead zone. The 
turbulence in the active zone propagates spiral density 
waves into the dead zone as seen in the snapshots at 
f /To = 5 (where To = 2k is the orbital period at ro). 

In Fig. 2 we show that the conditions for the devel- 
opment of the RWI are fulfilled in this disk model. We 
plot the azimuthal and vertical averages as a function of 
radius of density, azimuthal velocity, vorticity, and mag- 
netic field, for the first three orbits. In a barotropic disk, 
an extremum in potential vorticity cv z /p will launch the 
RWI (Lovelace et al. 1999). This condition is brought 
about by the transition in magnetic pressure near r n . 

At later times (lower panels of Fig. 1), a giant vortex 
is seen on the dead side of the transition. The density 



enhancement (lower left) matches spatially a vorticity 
minimum, confirming its anticyclonic nature. 

We show in Fig. 3 the time evolution of the azimuthal 
and vertical averages of the same quantities shown in 
Fig. 1. In the density plot we see the vortex as an en- 
hancement at the dead size of r n . It weakens with time 
but attains a steady state after £/To=60. The vorticity 
plot shows that the density enhancement is traced by 
anticyclonic vorticity. Conversely the density dips are 
spatially correlated with regions of cyclonic vorticity. 
Though the resistivity impedes the growth of the MRI 
in the dead zone, it experiences diffusion of some mag- 
netic field from the active zone. It thus appears weakly 
magnetized instead of completely demagnetized. 

Interestingly, in the same plots we see that vortices 
are excited in the magnetized regions. A weak vortex 
is seen in the middle of the active zone at very early 
times, t/To < 5. It eventually decays, after 20 orbits. 
Other intermittent vortices are also seen to be excited, 
at r=0.7, between 40 and 70 orbits, and again after 85 or- 
bits. These vortices exist in the midst of MRI turbulence, 
as in the models of Fromang & Nelson (2005). Strictly 
speaking, Lyra & Klahr (2011) state that vortices that are 
excited by the baroclinic instability and survive in hy- 
drodynamical models are destroyed when a magnetic 
field is abruptly introduced in the simulation. They do 
not exclude excitation of MHD vortices by other ways. 
Still, these magnetized vortices should host magneto- 
elliptic instabilities (Mizerski & Bajer 2009; Mizerski & 
Lyra 2012), and given their strength, should not exist; 
unless vorticity injection by the RWI is faster than de- 
struction by the magneto-elliptic instability in the rele- 
vant scales. 

The origin of the vortex at r=0.7 lies in the initial con- 
dition, because we used a field with a peak in magnetic 
pressure around that radius. This is seen in Fig. 2 as the 
peak in magnetic pressure translates into an azimuthal 
velocity inversion and a peak of potential vorticity. This 
incurs a non-equilibrium configuration, that induces 
sub-Keplerian motion inwards of the peak, and super- 
Keplerian motion outwards. It triggers a localized zonal 
flow, that in turns excites the non-axisymmetric modes 
of the RWI. 

3.2. Angular momentum transport 

We display in Fig. 4 the measured Maxwell and 
Reynolds stresses. The Reynolds stress is shown for 
both the active and dead zone. We see that the Maxwell 
and Reynolds stresses in the active zone reach the 10~ 2 
level, agreeing with previous studies (Papaloizou & 
Nelson 2003; Fromang & Nelson 2006; Lyra et al. 2008a; 
Dzyurkevich et al. 2010; Beckwith et al. 2011; Flock et al. 
2011; Sorathia et al. 2012). The novel result of our model 
is that the Reynolds stress in the dead zone is also at the 
10~ 2 level, only slightly lower than in the active zone. 

Without viscous dissipation, the angular momen- 
tum is transported to the outer boundary, where it is 
damped, while the matter accretes. As a result of out- 
ward angular momentum transport, the vortex should 
migrate inwards (Paardekooper et al. 2010; Lyra & Klahr 
2011) However, this is not observed in the simulations. 
The vortex keeps radiating waves, yet it sits in the pres- 
sure maximum that generated it, which acts as a migra- 
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FIG. 4. — Maxwell and Reynolds stresses as a function of time in the active and dead zones. Measurements were taken every To = 2n, and shown 
in the plot as box-and-whiskers five point summaries. The light thin grey lines are the whiskers marking the minimum and maximum values, the 
dark grey thick lines marks the lower and upper quartiles that box 50% of the values. A thicker colored line traces the median. High levels of 
Reynolds stress are maintained in the dead zone by the spiral density waves excited by the vortex at the active/ dead boundary. 
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FIG. 5.— Models of resolution 192x384, 384x768 and 768x1536 in 
the midplane. The vertical resolution is 64 in the first two plots, and 
128 in the last one. The vortex in the magnetized zone, that should 
be unstable, shrinks in size with increasing resolution, but is not 
quenched even in the highest resolution model (see text). The outer 
vortex exists in the dead zone and is not subject to the same destabi- 
lizing mechanism. The same colorbar is used as in Fig. 1. 



tion barrier. The same is seen in the models of Meheut 
et al. (2010, 2012b). 

3.3. Resolution study 

According to the results of Lyra & Klahr (2011), the 
magneto-elliptic instability should not allow the forma- 
tion of the vortices seen in the MRI-active region. This 
suggests that the models presented are under-resolving 
the MEI-modes. To examine this possibility we run 
models with twice (N r = 384, Afy = 768, N z = 64) and 
four times (N r = 768, N<p = 1536, N z = 128) the midplane 
resolution of the fiducial model. The flow pattern after 
12 orbits is shown in Fig. 5. We see that the inner vortex 
in the magnetized region has shrunk in size while the 
outer one in the dead zone has not. This is evidence that 
the MEI is being excited, and the vortex core shrinks to 
a size where it is under-resolved. The size of the scale 
height in all models is H = hr = 0.1 r, so H=0.08 for the 



location of the inner vortex. The low resolution model 
resolves that scale with H/ Ar=9 points, the middle res- 
olution one with 19 points and the highest resolution 
model with 38 points. 

We see that even in the highest resolution model, the 
density enhancement persists. We conclude that this 
vortex is either magnetically stable (which is problem- 
atic in view of the claims of Lyra & Klahr 2011) or the 
resolution requirement to capture the unstable modes 
has not been met. As we could not push the resolution 
even further (the highest resolution model already con- 
sumed 2 million computer hours), the question of the 
stability of the inner vortex will be addressed in a fu- 
ture work, in local box models, to satisfy the resolution 
requirement. 

3.4. 2tz model 

As shown by Flock et al. (2011) the main features 
of the MRI can be captured with a domain spanning 
7r in azimuth as well as they can in a domain span- 
ning 2/T. That statement, however, though applying for 
the axisymmetric MRI, does not automatically apply to 
the non-axisymmetric RWI, for which the most unstable 
wavelengths are m = 4 and m = 5 (Li et al. 2000). Vis- 
cous mergings then transfer the power towards lower 
wavenumbers, in a cascade leading to m = 1, i.e., a sin- 
gle vortex 6 . 

We check if there are significant differences pertain- 
ing to azimuthal domain size by performing a simu- 
lation spanning 2n radians and comparing the result 
with the fiducial model. The model has the same nu- 
merical resolution, using twice the grid points in <p, 
[N r ,N^,N z ]=[192,768,64]. The comparison is shown in 
Fig. 6. It is seen that at 12 orbits the In model shows a 
conspicuous m=2 mode, whereas the n model contains 
a single vortex at the same time. The two vortices in the 
2n model later merge into am = 1 mode, and as a result 

6 The m = 1 mode can cascade further, back to m = 0, thus turning 
the RWI into a cycle (see Regaly et al. 2012). However, this develop- 
ment occurs only over long evolution times and will not be dealt with 
in this work. 
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FIG. 6. — Since the RWI is non-axisymmetric, we assess the influence of azimuthal domain size Lq on its evolution by comparing models of 
Lq = n (left panels) and L^,=2n (right panels). Upper panels show the models at 12 orbits, lower panels at 26 orbits. Despite differences at early 
times (at 12 orbits the vortex cascade is still at m = 2 in the In models), once the inverse cascade is complete and a single vortex survives, the 
models look remarkably similar. 




FIG. 7. — Same as Fig. 4 but for the full In model. The level of angular momentum transport is very similar to the model spanning only n in 
azimuth. We conclude that the latter model leaves enough azimuthal room to avoid vortex self-interaction across the periodic boundary. 
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the flows in both cases are very similar at £=26 orbits. 

We also measure the angular momentum transport in 
the In model, shown in Fig. 7. The stresses are very sim- 
ilar to the ones shown in Fig. 4, for the n model. This 
means a box size of h^=n already provides enough az- 
imuthal spacing so that the vortex does not interact with 
itself in a way that significantly affects the stresses (the 
same is not true for local boxes). 

We conclude that studies concerning the linear 
growth of the instability require full In simulations, yet 
7r suffices if the interest is in the saturated, m — 1, state. 

4. CONCLUSIONS 

We performed the first simulation of a RWI- 
unstable boundary between dead and active 
zones in a protoplanetary disk with 3D resistive- 
magnetohydrodynamical models. We use a sharp 
resistive transition, that should be characteristic of the 
inner dead zone boundary. 

We find that Rossby vortices are excited on the dead 
side of the transition. They first appear in an m = 4 
mode, as expected from the RWI, then quickly merge 
viscously into a m = 2 mode, and finally reach m — 1, 
becoming a single giant vortex immediately outside of 
the transition. This lends credence to previous works 
that relied on alpha disk models with angular velocity 
transitions to trigger the RWI (Inaba & Barge 2006; Lyra 
et al. 2008b, 2009; Meheut et al. 2010, 2012b). 

We also assess the mass accretion rate in the dead 
zone, which is a result both of waves propagating from 
the active zone and of waves radiated by the vortex. We 
find normalized Reynolds stresses of similar magnitude 
in both the active and in the dead zones, at the 10~ 2 
level. Vortex-free transitions develop dead zone stresses 
two orders of magnitude quieter than the active zones, 
so this high accretion rate is a manifestation of the abil- 
ity of the vortex to maintain high accretion rates. The 
same level is found for angular momentum transport in 
baroclinically unstable local disk models (Lesur & Pa- 
paloizou 2009; Lyra & Klahr 2011). Because the vortex is 
born in a high pressure ridge, it does not migrate further 
in; the ridge acts as a migration barrier. 

A puzzling aspect of our model is the development 
of a weak vortex in the middle of the active zone. The 
vortex was driven by the RWI at the peak in magnetic 
pressure due to the artificial initial condition. Yet ac- 
cording to a recent study by Lyra & Klahr (2011), ro- 
tating magnetized vortices should be host to a set of 
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unstable centrifugal hydromagnetic modes (magneto- 
elliptic instability) that should ultimately break closed 
elliptic streamlines. If that is the case, then our mod- 
els are either under-resolving these modes, or vorticity 
injection by the RWI is faster than destruction by the 
magneto-elliptic instability. We pushed the resolution 
of the global model to the limit allowed by current com- 
putational power. The vortex shrank with increasing 
resolution, but we could not quench it completely. A 
conclusive test will require local models, which we plan 
to undertake in future work. A study of convergence 
with dead zone size, and an assessment of how wide 
the resistive jump can be while still allowing for Rossby 
vortex excitation are also lacking. 

The models presented leave room for exciting de- 
velopments. First, we used a static resistivity profile. 
Though this is shown by Latter & Balbus (2012) to pro- 
duce acceptable results, we plan to include dynamic re- 
sistivity in future models, following the reduced recom- 
bination network of Ilgner & Nelson (2006a,b,c), used 
by Turner & Sano (2008); Turner & Drake (2009); Gres- 
sel et al. (2011); Hirose & Turner (2011) and Gressel 
et al. (2012). Effects of Hall MHD and ambipolar dif- 
fusion (Wardle 1999; Salmeron & Wardle 2005; Bai & 
Stone 2011; Bai 2011) may also impact significantly on 
the nature of the active / dead radial transition. Future 
work will also include stratification, shown through nu- 
merical (Meheut et al. 2010, 2012b) and analytical work 
(Meheut et al. 2012a; Lin 2012) to matter significantly 
for the internal dynamics of 3D vortices. Finally, inter- 
acting particles should be included to assess the ability 
of three-dimensional Rossby vortices in dead zones to 
assist on planet formation. 
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